## load packages
library(mgcv)
library(plyr)
library(visreg)
library(ggplot2)
library(MuMIn)
library(reticulate)

Model Building

## read in the main data
allom_data <- read.csv('~/Desktop/References_allometry/WMNF_allom5.csv', header = TRUE)
str(allom_data)
## 'data.frame':    318 obs. of  21 variables:
##  $ SPP4       : Factor w/ 8 levels "ACPE","ACRU",..: 2 2 2 3 3 3 4 4 4 5 ...
##  $ ID         : Factor w/ 68 levels "1","10","11",..: 55 55 55 55 55 55 55 55 55 55 ...
##  $ ELEVB      : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ TREEAGE    : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ STANDAGE   : int  14 14 14 14 14 14 14 14 14 14 ...
##  $ LOCATION   : Factor w/ 3 levels "BEF","HBEF","WMNF": 1 1 1 1 1 1 1 1 1 1 ...
##  $ DataSetID  : Factor w/ 6 levels "BattlesGB","BattlesRock",..: 4 4 4 4 4 4 4 4 4 4 ...
##  $ SOURCE     : Factor w/ 5 levels "BattlesJ","FaheyT",..: 3 3 3 3 3 3 3 3 3 3 ...
##  $ DBHcm      : num  2.2 3.4 4.4 2 2.9 3.9 2.4 3 3.8 3.3 ...
##  $ HTcm       : num  554 520 685 466 647 569 455 430 625 840 ...
##  $ PVcm3      : num  1053 2361 5208 732 2137 ...
##  $ LEAFg      : int  69 179 415 89 279 341 97 280 330 142 ...
##  $ TWIGg      : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ CTLGRg     : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ LIVEBRANCHg: num  281.5 376.8 779.9 97.4 311.6 ...
##  $ BOLEBARKg  : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ BOLEWOODg  : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ BOLETOTg   : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ ABOVEg     : num  848 1986 4283 745 1916 ...
##  $ DEADWOODg  : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ TOTROOTg   : num  NA NA NA NA NA NA NA NA NA NA ...
allom_data[, "STANDAGE"] <- as.factor(allom_data[, "STANDAGE"])
## Remove the data from the 98 year old stand so that all of the data are comparable:
allom_data_Y <- subset(allom_data, STANDAGE != "98years")
str(allom_data_Y)
## 'data.frame':    318 obs. of  21 variables:
##  $ SPP4       : Factor w/ 8 levels "ACPE","ACRU",..: 2 2 2 3 3 3 4 4 4 5 ...
##  $ ID         : Factor w/ 68 levels "1","10","11",..: 55 55 55 55 55 55 55 55 55 55 ...
##  $ ELEVB      : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ TREEAGE    : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ STANDAGE   : Factor w/ 7 levels "14","16","17",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ LOCATION   : Factor w/ 3 levels "BEF","HBEF","WMNF": 1 1 1 1 1 1 1 1 1 1 ...
##  $ DataSetID  : Factor w/ 6 levels "BattlesGB","BattlesRock",..: 4 4 4 4 4 4 4 4 4 4 ...
##  $ SOURCE     : Factor w/ 5 levels "BattlesJ","FaheyT",..: 3 3 3 3 3 3 3 3 3 3 ...
##  $ DBHcm      : num  2.2 3.4 4.4 2 2.9 3.9 2.4 3 3.8 3.3 ...
##  $ HTcm       : num  554 520 685 466 647 569 455 430 625 840 ...
##  $ PVcm3      : num  1053 2361 5208 732 2137 ...
##  $ LEAFg      : int  69 179 415 89 279 341 97 280 330 142 ...
##  $ TWIGg      : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ CTLGRg     : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ LIVEBRANCHg: num  281.5 376.8 779.9 97.4 311.6 ...
##  $ BOLEBARKg  : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ BOLEWOODg  : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ BOLETOTg   : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ ABOVEg     : num  848 1986 4283 745 1916 ...
##  $ DEADWOODg  : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ TOTROOTg   : num  NA NA NA NA NA NA NA NA NA NA ...
max(allom_data_Y$DBHcm)
## [1] 66
## get subset of data which includes trees in 98year stand which are smaller than the ## largest trees in the younger stands.

allom_data_trunc <- subset(allom_data, DBHcm <= 66)
## log transform aboveground biomass:
allom_data_trunc$logABOVEg <- log(allom_data_trunc$ABOVEg)
ggplot(allom_data_trunc, aes(x = DBHcm, y = logABOVEg)) +
  geom_jitter() +
  geom_smooth(size = 0.5)

## lets try taking the log transformation of both variables:
allom_data_trunc$logDBHcm <- log(allom_data_trunc$DBHcm)
ggplot(allom_data_trunc, aes(x = logDBHcm, y = logABOVEg)) +
  geom_jitter() +
  geom_smooth(size = 0.5)

## Create SPPSTAGE in order to fully cross in GAM
for (i in 1:length(allom_data_trunc$SPP4)) {
  allom_data_trunc$SPPSTAGE[i] <- paste(allom_data_trunc$SPP4[i], allom_data_trunc$STANDAGE[i])
}

allom_data_trunc$SPPSTAGE <- as.factor(allom_data_trunc$SPPSTAGE)


## log transform HTcm:
allom_data_trunc$logHTcm <- log(allom_data_trunc$HTcm)

## log transform PVcm3
allom_data_trunc$logPVcm3 <- log(allom_data_trunc$PVcm3)

## remove NA values of biomass:
allom_data_trunc <- subset(allom_data_trunc, logABOVEg != "NA")
## look at how many data points we have for each species.  ACRU and PIRU seem low. 
## Let's create a dataset with them removed to compare to one with them included.

count(allom_data_trunc$SPP4)
##      x freq
## 1 ACPE   40
## 2 ACRU   12
## 3 ACSA   49
## 4 BEAL   62
## 5 BEPA   41
## 6 FAGR   48
## 7 PIRU   15
## 8 PRPE   49
allom_data_truncSPP4 <- subset(allom_data_trunc, SPP4 != "ACRU")
allom_data_truncSPP4 <- subset(allom_data_truncSPP4, SPP4 != "PIRU")
count(allom_data_truncSPP4$SPP4)
##      x freq
## 1 ACPE   40
## 2 ACSA   49
## 3 BEAL   62
## 4 BEPA   41
## 5 FAGR   48
## 6 PRPE   49
## make another subset of the data with only SPPSTAGE levels that have more than 10 data points.

SPPSTAGE_counts <- data.frame(count(allom_data_trunc$SPPSTAGE))
SPPSTAGE_counts <- subset(SPPSTAGE_counts, freq >= 10)

SPPSTAGE_list <- as.factor(SPPSTAGE_counts$x)

allom_data_truncSPPSTAGE <- subset(allom_data_trunc, (allom_data_trunc$SPPSTAGE %in% SPPSTAGE_list) == TRUE)

allom_data_truncSPPSTAGE$SPPSTAGE
##   [1] BEAL 17 BEAL 17 BEAL 17 BEAL 17 BEAL 17 BEAL 17 BEPA 17 BEPA 17
##   [9] BEPA 17 BEPA 17 BEPA 17 BEPA 17 PRPE 17 PRPE 17 PRPE 17 PRPE 17
##  [17] PRPE 17 PRPE 17 PRPE 17 PRPE 17 PRPE 17 PRPE 17 PRPE 17 PRPE 17
##  [25] PRPE 17 PRPE 17 PRPE 17 ACPE 17 ACPE 17 ACPE 17 ACPE 17 ACPE 17
##  [33] ACPE 17 ACPE 17 ACPE 17 ACPE 17 ACPE 17 ACSA 17 ACSA 17 ACSA 17
##  [41] ACSA 17 ACSA 17 ACSA 17 ACSA 17 ACSA 17 ACSA 17 ACSA 17 ACSA 17
##  [49] BEAL 17 BEAL 17 BEAL 17 BEAL 17 BEAL 17 BEAL 17 BEAL 17 BEAL 17
##  [57] BEAL 17 BEAL 17 BEAL 17 BEPA 17 BEPA 17 BEPA 17 BEPA 17 BEPA 17
##  [65] BEPA 17 BEPA 17 BEPA 17 BEPA 17 BEPA 17 BEPA 17 FAGR 17 FAGR 17
##  [73] FAGR 17 FAGR 17 FAGR 17 FAGR 17 FAGR 17 FAGR 17 FAGR 17 FAGR 17
##  [81] FAGR 17 PRPE 17 PRPE 17 PRPE 17 PRPE 17 PRPE 17 PRPE 17 PRPE 17
##  [89] PRPE 17 PRPE 17 PRPE 17 ACPE 26 ACPE 26 ACPE 26 ACPE 26 ACPE 26
##  [97] ACPE 26 ACPE 26 ACPE 26 ACPE 26 ACPE 26 ACPE 26 ACPE 26 ACPE 26
## [105] ACPE 26 ACPE 26 BEAL 26 BEAL 26 BEAL 26 BEAL 26 BEAL 26 BEAL 26
## [113] BEAL 26 BEAL 26 BEPA 26 BEPA 26 BEPA 26 BEPA 26 BEPA 26 BEPA 26
## [121] BEPA 26 BEPA 26 PRPE 26 PRPE 26 PRPE 26 PRPE 26 PRPE 26 PRPE 26
## [129] PRPE 26 PRPE 26 BEAL 26 BEAL 26 BEAL 26 BEPA 26 BEPA 26 BEPA 26
## [137] PRPE 26 PRPE 26 PRPE 26 ACPE 98 ACPE 98 ACPE 98 ACPE 98 ACPE 98
## [145] ACPE 98 ACPE 98 ACPE 98 ACPE 98 ACPE 98 ACPE 98 ACPE 98 ACPE 98
## [153] ACPE 98 ACPE 98 ACSA 98 ACSA 98 ACSA 98 ACSA 98 ACSA 98 ACSA 98
## [161] ACSA 98 ACSA 98 ACSA 98 ACSA 98 ACSA 98 ACSA 98 ACSA 98 ACSA 98
## [169] ACSA 98 ACSA 98 ACSA 98 ACSA 98 ACSA 98 ACSA 98 ACSA 98 BEAL 98
## [177] BEAL 98 BEAL 98 BEAL 98 BEAL 98 BEAL 98 BEAL 98 BEAL 98 BEAL 98
## [185] BEAL 98 BEAL 98 BEAL 98 BEAL 98 BEAL 98 BEAL 98 BEAL 98 BEAL 98
## [193] BEAL 98 BEAL 98 BEAL 98 BEAL 98 FAGR 98 FAGR 98 FAGR 98 FAGR 98
## [201] FAGR 98 FAGR 98 FAGR 98 FAGR 98 FAGR 98 FAGR 98 FAGR 98 FAGR 98
## [209] FAGR 98 FAGR 98 FAGR 98 FAGR 98 FAGR 98 FAGR 98 FAGR 98 FAGR 98
## [217] FAGR 98 PIRU 98 PIRU 98 PIRU 98 PIRU 98 PIRU 98 PIRU 98 PIRU 98
## [225] PIRU 98 PIRU 98 PIRU 98 PIRU 98 PIRU 98 PIRU 98 PIRU 98 PIRU 98
## 39 Levels: ACPE 17 ACPE 26 ACPE 98 ACRU 14 ACRU 16 ACRU 26 ... PRPE 29
## Try STANDAGE as numerical

allom_data_truncSTANDAGENUM <- allom_data_trunc
allom_data_truncSTANDAGENUM$STANDAGE <- as.numeric(allom_data_trunc$STANDAGE)

Perhaps in the future we could keep the small trees only from the 98 year data as there is some basis to comparison, although its perhaps still too large of a gap to justify any conclusions.

It looks like we get something roughly linear with the loglog transformation of DBH vs BIOMASS. Based on a visual inspection alone, it appears that the heteroscedasticity issue has been at least significantly reduced.

Define Function to Create List of Models (linear and GAM):

A note: A more experienced/clever programmer could probably do this in 1/4 the lines of code I used. I am no such programmer, so this is likely the most overcomplicated and inefficient method of doing this. Oh well, maybe some day…

## comparing all possible combinations of linear models for the dataset using AIC and 
## AICc:
## knowledge of the dredge function may have been useful/done the work for me for the gam 
## models
options(na.action = "na.fail")
model1 <- lm(logABOVEg ~ logDBHcm * logHTcm * SPP4 * STANDAGE, data = allom_data_trunc)
model_sel_linear_AIC <- dredge(model1, rank = AIC)
model_sel_linear_AICc <- dredge(model1, rank = AICc)


## Function to write out all possible model terms, to be combined in later funcitons:
build_model_terms <- function(x, by, op) {
  s_terms <- character()
  int_terms <- character()
  ## for loops to write all possible smooth terms:
  ## l is dummy variable specifying inclusion of id = 1:
  for (l in 1:2) {
    if (l == 1) {
      ## k is a dummy variable indicating whether or not theres an interaction term:
      for (k in 1:2) {
        if (k == 1) {
          for (i in 1:length(x)) {
            for (j in 1:length(by)) {
              s_termtext <- paste0("s(", x[i], ", k = 4", ", by =", by[j], " ,id = 1", ")")
              s_terms <- c(s_terms, s_termtext)
              }
            } 
          }
        else if (k == 2) {
          for (j in 1:length(by)) {
            s_termtext <- paste0("s(", x[1], ", ", x[2], ", k = 4", ", by =", by[j], " ,id = 1", ")")
            s_terms <- c(s_terms, s_termtext)
            }
          }
        }
      }
    else if (l == 2) {
      ## k is a dummy variable indicating whether or not theres an interaction term:
      for (k in 1:2) {
        if (k == 1) {
          for (i in 1:length(x)) {
            for (j in 1:length(by)) {
              s_termtext <- paste0("s(", x[i], ", k = 4", ", by =", by[j], ")")
              s_terms <- c(s_terms, s_termtext)
              }
            } 
          }
        else if (k == 2) {
          for (j in 1:length(by)) {
            s_termtext <- paste0("s(", x[1], ", ", x[2], ", k = 4", ", by =", by[j], ")")
            s_terms <- c(s_terms, s_termtext)
            }
          }
        }
      }
    }
  ## k is a dummy variable indicating whether there are 1 or 2 intercept terms
  for (k in 1:2) {
      if (k == 1) {
        for (i in 1:length(by)) {
          by1 <- by[! by %in% " "]
          if (!is.na(by1[i])) {
            int_termtext <- paste0(by[i])
            int_terms <- c(int_terms, int_termtext)
          }
          else {}
        }  
      }
      else if (k == 2) {
        for (l in 1:2) {
          int_termtext <- paste0(by[1], op[l], by[2])
          int_terms <- c(int_terms, int_termtext)
        }
      }
    }
  model_terms <- list(int_terms, s_terms)
  return(model_terms)
  }


## Function builds models with one s() term - including those with two explanatory variables and an 
## interaction:
build_models_1 <- function(y, x, by, int_terms, s_terms) {
  models <- character()
  for (i in 1:length(s_terms)) {
    for (j in 1:length(int_terms)) {
      ## k is a dummy variable that indicates whether there are 1 or two intercept terms
      for (k in 1:2) {
        if (k == 1) {
          ## When we have only 1 intercept term:
          if (!grepl("\\+", int_terms[j]) && !grepl("\\*", int_terms[j])) {
            if (grepl(int_terms[j], s_terms[i]) && !grepl("SPPSTAGE", s_terms[i])) {
              newmodeltext <- paste0(y, "~ ", s_terms[i], " + ", int_terms[j])
              newmodel <- as.formula(parse(text = newmodeltext)[[1]])
              models <- c(models, newmodel)
            }
            else if (grepl("SPPSTAGE", s_terms[i])) {
              newmodeltext <- paste0(y, " ~ ", s_terms[i], " + ", "SPPSTAGE")
              newmodel <- as.formula(parse(text = newmodeltext)[[1]])
              models <- c(models, newmodel)
            }
            else if (grepl("by = )", s_terms[i]) && !grepl("SPPSTAGE", int_terms[j])) {
              newmodeltext <- paste0(y, " ~ ", s_terms[i], " + ", int_terms[j])
              newmodel <- as.formula(parse(text = newmodeltext)[[1]])
              models <- c(models, newmodel)
            }
            else {}
          }  
          else {}
        }
        else if (k == 2) {
          if (grepl("SPPSTAGE", s_terms[i])) {}
          else if (grepl("\\+", int_terms[j]) | grepl("\\*", int_terms[j])) {
            newmodeltext <- paste0(y, " ~ ", s_terms[i], " + ", int_terms[j])
            newmodel <- as.formula(parse(text = newmodeltext)[[1]])
            models <- c(models, newmodel)
          }
        }
      }
    }  
  }
  models <- unique(models)
  return(models)
}

## Function to build models for case with two s() terms, no interactions:
## I got a bit carried away with this one I think.
build_models_2 <- function(y, x, by, int_terms, s_terms) {
  models <- character()
  ## create a vector of DBH s() terms and one of HT s() terms
  dbh_terms <- character()
  ht_terms <- character()
  for (i in 1:length(s_terms)) {
    if (grepl("logDBHcm", s_terms[i]) && !grepl("logHTcm", s_terms[i])) {
      dbh_terms <- c(dbh_terms, s_terms[i])
    }
    else if (grepl("logHTcm", s_terms[i]) && !grepl("logDBHcm", s_terms[i])) {
      ht_terms <- c(ht_terms, s_terms[i])
    }
    else {}
  }
  ## Make sure there are no duplicates
  dbh_terms <- unique(dbh_terms)
  ht_terms <- unique(ht_terms)
  
  for (i in 1:length(dbh_terms)) {
    for (l in 1:length(ht_terms)) {
      for (j in 1:length(int_terms)) {
        for (k in 1:2) {
          if (k == 1) {
            ## When we have only 1 intercept term (except with special case of SPP4STANDAGE):
            if (!grepl("\\+", int_terms[j]) && !grepl("\\*", int_terms[j])) {
              ## for the case without SPP4STANDAGE, require that both dbh_term and ht_term are represented in int_term:
              if ((grepl(int_terms[j], dbh_terms[i]) && grepl(int_terms[j], ht_terms[l])) | (grepl(int_terms[j], paste0(dbh_terms[i], ht_terms[l])) && grepl("by = )", paste0(dbh_terms[i], ht_terms[l]))) && !grepl("SPPSTAGE", paste0(dbh_terms[i], ht_terms[l]))) {
                newmodeltext <- paste0(y, " ~ ", dbh_terms[i], " + ", ht_terms[l], " + ", int_terms[j])
                newmodel <- as.formula(parse(text = newmodeltext)[[1]])
                models <- c(models, newmodel)
              }
              ## for the case with SPPSTAGE:
              else if (grepl("SPPSTAGE", paste0(dbh_terms[i], ht_terms[i]))) {
                ## for the case with SPP4STANDAGE and different by term, require that diff by term is represented in int_term
                if (grepl(int_terms[j], paste0(dbh_terms[i], ht_terms[l])) && int_terms[j] != "SPPSTAGE") {
                  newmodeltext <- paste0(y, " ~ ", dbh_terms[i], " + ", ht_terms[l], " + ", "SPPSTAGE", " + ", int_terms[j])
                  newmodel <- as.formula(parse(text = newmodeltext)[[1]])
                  models <- c(models, newmodel)
                }
                ## for the case with only SPPSTAGE as a by term, only SPPSTAGE as an int_term
                else if (grepl(int_terms[j], dbh_terms[i]) && grepl(int_terms[j], ht_terms[l])) {
                  newmodeltext <- paste0(y, " ~ ", dbh_terms[i], " + ", ht_terms[l], " + ", "SPPSTAGE")
                  newmodel <- as.formula(parse(text = newmodeltext)[[1]])
                  models <- c(models, newmodel)
                }
                else {}
              }
              ## for the case with empty by terms
              else if (grepl("by = )", dbh_terms[i]) && grepl("by = )", ht_terms[l]) && !grepl("SPPSTAGE", int_terms[j])) {
                newmodeltext <- paste0(y, " ~ ", dbh_terms[i], " + ", ht_terms[l], " + ", int_terms[j])
                newmodel <- as.formula(parse(text = newmodeltext)[[1]])
                models <- c(models, newmodel)
              }
              else {}
            }  
            else {}
          }
          else if (k == 2) {
            if (grepl("SPPSTAGE", paste0(int_terms[j], dbh_terms[i], ht_terms[l]))) {}
            else if (grepl("\\+", int_terms[j]) | grepl("\\*", int_terms[j])) {
              newmodeltext <- paste0(y, " ~ ", dbh_terms[i], " + ", ht_terms[l], " + ", int_terms[j])
              newmodel <- as.formula(parse(text = newmodeltext)[[1]])
              models <- c(models, newmodel)
            }
          }
        }
      }  
    }
  }
  models <- unique(models)
  return(models)
}

## Function to put it all together:
allom_gams <- function(data1, y, x, by, rank = "AIC") {
  ## y, x and by must be input as concatenated string variables
  x <- c("logDBHcm", "logHTcm", "logPVcm3")
  by <- c("SPP4", "STANDAGE", "SPPSTAGE")
  by <- append(by, " ")
  op <- c(" + ", " * ")
  
  ## call function to build models for case with one explanatory variable
  model_terms <- build_model_terms(x = x, by = by, op = op)
  int_terms <- model_terms[[1]]
  s_terms <- model_terms[[2]]
  models_1 <- build_models_1(y = y, x = x, by = by, int_terms = int_terms, s_terms = s_terms)
  models_2 <- build_models_2(y = y, x = x, by = by, int_terms = int_terms, s_terms = s_terms)
  models <- c(models_1, models_2)
  
  ## fit gam models:
  for (i in 1:length(models)) {
    assign(paste0("model", i), gam(models[[i]], data = data1))
  }
  
  model_list <- character()
  for (i in 1:length(models)) {
    if (i < max(length(models))) {
      model_list <- paste0(model_list, "model", i, ", ")
    }
    else {
      model_list <- paste0(model_list, "model", i)
    }
  }
  results <- eval(parse(text = paste("model.sel(", model_list, ", rank =", rank, ")")))
  return(results)   
}

## build, run and perform model selection on models
model_sel_gam_AIC <- allom_gams(data1 = allom_data_trunc, y = "logABOVEg", x = c("logDBHcm", "logHTcm", "logPVcm3"), by = c("SPP4", "STANDAGE", "SPPSTAGE"), rank = "AIC")

## Check model selection on dataset without poorly represented species to make sure results don't change:
# model_sel_gam_AICtruncSPP4 <- allom_gams(data1 = allom_data_truncSPP4, y = "logABOVEg", x = c("logDBHcm", "logHTcm", "logPVcm3"), by = c("SPP4", "STANDAGE", "SPPSTAGE"), rank = "AIC")
# 
# ## Check model selection on dataset without poorly represented species/standage combos:
# model_sel_gam_AICtruncSPPSTAGE <- allom_gams(data1 = allom_data_truncSPPSTAGE, y = "logABOVEg", x = c("logDBHcm", "logHTcm", "logPVcm3"), by = c("SPP4", "STANDAGE", "SPPSTAGE"), rank = "AIC")
# 
# ## Check model selection when STANDAGE treated as numeric instead of factor:
# model_sel_gam_AICtruncSTANDAGENUM <- allom_gams(data1 = allom_data_truncSTANDAGENUM, y = "logABOVEg", x = c("logDBHcm", "logHTcm", "logPVcm3"), by = c("SPP4", "STANDAGE", "SPPSTAGE"), rank = "AIC")

Here are the results when using AIC. The output becomes quite cumbersome when linear and gam models are combined because of the number of possible terms between the two model families. The four models with the lowest AIC score are gam models, the fifth being linear.

Here they are listed out in a more readable format:

1) logABOVEg ~ s(logDBHcm, by = SPPSTAGE) + s(logHTcm, by = SPPSTAGE) + SPPSTAGE
2) logABOVEg ~ s(logDBHcm, by = SPPSTAGE) + s(logHTcm, by = SPP4) + SPPSTAGE
3) logABOVEg ~ s(logDBHcm, by = SPPSTAGE) + s(logHTcm, by = STANDAGE) + SPPSTAGE
4) logABOVEg ~ s(logPVcm3, by = SPPSTAGE) + SPPSTAGE
5) logABOVEg ~ logDBHcm + logHTcm + SPP4 + STANDAGE + logDBHcm:SPP4 + logDBH:STANDAGE +
logHT:SPP4 + logHT:STANDAGE + SPP4:STANDAGE + logHT:SPP4:STANDAGE
model_sel_results_AIC <- merge(model_sel_linear_AIC, model_sel_gam_AIC)
head(model_sel_results_AIC)
## Model selection table 
##            (Int) SPP4 STA s(lDB,4,SPPS,1) SPPS s(lHT,4,SPP4,1)
## model164.y 7.949                             +                
## model162.y 7.867    +                        +                
## model159.y 8.130    +                        +               +
## model116.y 7.949                        +    +                
## model163.y 7.961        +                    +                
## model161.y 7.969                             +                
##            s(lHT,4,SPPS,1) s(lDB,4,SPPS) s(lHT,4,SPP4) s(lHT,4,STA)
## model164.y                             +                           
## model162.y                             +             +             
## model159.y                             +                           
## model116.y                                                         
## model163.y                             +                          +
## model161.y               +             +                           
##            s(lHT,4,SPPS)  df  logLik   AIC delta weight
## model164.y             + 125 145.774 -41.0  0.00  0.681
## model162.y               103 122.893 -38.2  2.85  0.163
## model159.y               107 126.390 -37.4  3.62  0.111
## model116.y             + 123 141.314 -35.1  5.94  0.035
## model163.y                93 109.201 -32.0  9.03  0.007
## model161.y               129 144.594 -29.4 11.61  0.002
## Models ranked by AIC(x)
allom_model <- gam(logABOVEg ~ s(logDBHcm, by = SPPSTAGE, k = 4) + s(logHTcm, by = SPPSTAGE, k = 4) + SPP4 * STANDAGE, data = allom_data_trunc, qr=FALSE)

summary(allom_model)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## logABOVEg ~ s(logDBHcm, by = SPPSTAGE, k = 4) + s(logHTcm, by = SPPSTAGE, 
##     k = 4) + SPP4 * STANDAGE
## 
## Parametric coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           7.91407    0.55737  14.199   <2e-16 ***
## SPP4ACRU              1.06915    0.64682   1.653   0.1000 .  
## SPP4ACSA              1.28603    0.63987   2.010   0.0459 *  
## SPP4BEAL              1.20918    0.65695   1.841   0.0672 .  
## SPP4BEPA              1.22375    0.99314   1.232   0.2194    
## SPP4FAGR              1.71842    0.81252   2.115   0.0357 *  
## SPP4PIRU             -0.02006    0.13693  -0.146   0.8837    
## SPP4PRPE              0.90884    0.67657   1.343   0.1808    
## STANDAGE16           -5.71597    2.22182  -2.573   0.0108 *  
## STANDAGE17            0.03533    0.61209   0.058   0.9540    
## STANDAGE23           -1.61731    5.67910  -0.285   0.7761    
## STANDAGE26            0.96506    0.55979   1.724   0.0863 .  
## STANDAGE29            0.97583    4.09147   0.239   0.8117    
## STANDAGE98            0.97960    0.56124   1.745   0.0825 .  
## SPP4ACRU:STANDAGE16   5.60162    2.25370   2.486   0.0138 *  
## SPP4ACSA:STANDAGE16   4.66599    2.05804   2.267   0.0245 *  
## SPP4BEAL:STANDAGE16   4.25476    5.54187   0.768   0.4436    
## SPP4BEPA:STANDAGE16   4.61852    2.49501   1.851   0.0657 .  
## SPP4FAGR:STANDAGE16 -29.79298   12.04542  -2.473   0.0143 *  
## SPP4PIRU:STANDAGE16   0.00000    0.00000      NA       NA    
## SPP4PRPE:STANDAGE16   4.93611    2.26701   2.177   0.0307 *  
## SPP4ACRU:STANDAGE17   0.00000    0.00000      NA       NA    
## SPP4ACSA:STANDAGE17   4.35475    3.10772   1.401   0.1628    
## SPP4BEAL:STANDAGE17  -0.86693    0.78604  -1.103   0.2714    
## SPP4BEPA:STANDAGE17  -0.05194    1.02853  -0.050   0.9598    
## SPP4FAGR:STANDAGE17   8.50206   11.31072   0.752   0.4532    
## SPP4PIRU:STANDAGE17   0.00000    0.00000      NA       NA    
## SPP4PRPE:STANDAGE17   0.03935    0.72418   0.054   0.9567    
## SPP4ACRU:STANDAGE23   0.00000    0.00000      NA       NA    
## SPP4ACSA:STANDAGE23   0.00000    0.00000      NA       NA    
## SPP4BEAL:STANDAGE23   2.25911    5.78284   0.391   0.6965    
## SPP4BEPA:STANDAGE23   1.49034    5.73253   0.260   0.7952    
## SPP4FAGR:STANDAGE23   0.00000    0.00000      NA       NA    
## SPP4PIRU:STANDAGE23   0.00000    0.00000      NA       NA    
## SPP4PRPE:STANDAGE23  -5.36675   16.98443  -0.316   0.7524    
## SPP4ACRU:STANDAGE26  -1.34815    0.99896  -1.350   0.1788    
## SPP4ACSA:STANDAGE26  -1.01868    0.69343  -1.469   0.1435    
## SPP4BEAL:STANDAGE26  -0.93070    0.69367  -1.342   0.1813    
## SPP4BEPA:STANDAGE26  -1.31017    1.01748  -1.288   0.1994    
## SPP4FAGR:STANDAGE26  -1.33278    0.82017  -1.625   0.1058    
## SPP4PIRU:STANDAGE26   0.00000    0.00000      NA       NA    
## SPP4PRPE:STANDAGE26  -0.82925    0.73802  -1.124   0.2626    
## SPP4ACRU:STANDAGE29  -0.95117    4.11280  -0.231   0.8174    
## SPP4ACSA:STANDAGE29  -0.24322    4.20843  -0.058   0.9540    
## SPP4BEAL:STANDAGE29  -0.68423    4.11907  -0.166   0.8682    
## SPP4BEPA:STANDAGE29   4.95649   24.49107   0.202   0.8398    
## SPP4FAGR:STANDAGE29  -1.45498    4.12750  -0.353   0.7248    
## SPP4PIRU:STANDAGE29   0.00000    0.00000      NA       NA    
## SPP4PRPE:STANDAGE29  -0.64707    4.15499  -0.156   0.8764    
## SPP4ACRU:STANDAGE98   0.00000    0.00000      NA       NA    
## SPP4ACSA:STANDAGE98  -1.15923    0.64601  -1.794   0.0743 .  
## SPP4BEAL:STANDAGE98  -1.07396    0.66310  -1.620   0.1070    
## SPP4BEPA:STANDAGE98   0.00000    0.00000      NA       NA    
## SPP4FAGR:STANDAGE98  -1.49312    0.81746  -1.827   0.0693 .  
## SPP4PIRU:STANDAGE98  -0.02006    0.13693  -0.146   0.8837    
## SPP4PRPE:STANDAGE98   0.00000    0.00000      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                edf Ref.df       F  p-value    
## s(logDBHcm):SPPSTAGEACPE 17 1.0000 1.0000  12.602 0.000482 ***
## s(logDBHcm):SPPSTAGEACPE 26 1.0000 1.0000  28.241 2.81e-07 ***
## s(logDBHcm):SPPSTAGEACPE 98 1.0000 1.0000  28.734 2.25e-07 ***
## s(logDBHcm):SPPSTAGEACRU 14 1.0000 1.0000  14.472 0.000190 ***
## s(logDBHcm):SPPSTAGEACRU 16 1.0000 1.0000   6.096 0.014402 *  
## s(logDBHcm):SPPSTAGEACRU 26 1.0000 1.0000   2.735 0.099759 .  
## s(logDBHcm):SPPSTAGEACRU 29 1.0000 1.0000   4.961 0.027058 *  
## s(logDBHcm):SPPSTAGEACSA 14 1.0000 1.0000  18.927 2.17e-05 ***
## s(logDBHcm):SPPSTAGEACSA 16 0.4606 0.4606   0.025 0.915474    
## s(logDBHcm):SPPSTAGEACSA 17 1.0000 1.0000  10.772 0.001219 ** 
## s(logDBHcm):SPPSTAGEACSA 26 1.0000 1.0000  17.680 3.95e-05 ***
## s(logDBHcm):SPPSTAGEACSA 29 1.0000 1.0000   1.589 0.208941    
## s(logDBHcm):SPPSTAGEACSA 98 1.4033 1.6758 155.808  < 2e-16 ***
## s(logDBHcm):SPPSTAGEBEAL 14 1.0000 1.0000   8.167 0.004727 ** 
## s(logDBHcm):SPPSTAGEBEAL 16 1.0000 1.0000   0.038 0.845011    
## s(logDBHcm):SPPSTAGEBEAL 17 1.0000 1.0000  14.943 0.000150 ***
## s(logDBHcm):SPPSTAGEBEAL 23 1.0000 1.0000   1.524 0.218468    
## s(logDBHcm):SPPSTAGEBEAL 26 1.0000 1.0000  13.976 0.000242 ***
## s(logDBHcm):SPPSTAGEBEAL 29 1.0000 1.0000   4.091 0.044488 *  
## s(logDBHcm):SPPSTAGEBEAL 98 1.4009 1.6002 197.387  < 2e-16 ***
## s(logDBHcm):SPPSTAGEBEPA 14 1.0000 1.0000   2.361 0.126018    
## s(logDBHcm):SPPSTAGEBEPA 16 1.0000 1.0000   0.255 0.614217    
## s(logDBHcm):SPPSTAGEBEPA 17 1.0000 1.0000  50.511 1.67e-11 ***
## s(logDBHcm):SPPSTAGEBEPA 23 1.0000 1.0000  26.291 6.86e-07 ***
## s(logDBHcm):SPPSTAGEBEPA 26 1.0000 1.0000  51.077 1.31e-11 ***
## s(logDBHcm):SPPSTAGEBEPA 29 1.0000 1.0000   0.078 0.780913    
## s(logDBHcm):SPPSTAGEFAGR 14 1.0000 1.0000   6.648 0.010661 *  
## s(logDBHcm):SPPSTAGEFAGR 16 1.0000 1.0000   5.996 0.015222 *  
## s(logDBHcm):SPPSTAGEFAGR 17 1.6803 1.8858  23.388 2.36e-07 ***
## s(logDBHcm):SPPSTAGEFAGR 26 1.0000 1.0000  39.876 1.59e-09 ***
## s(logDBHcm):SPPSTAGEFAGR 29 0.3243 0.3243  61.359 1.36e-05 ***
## s(logDBHcm):SPPSTAGEFAGR 98 1.0000 1.0000 361.507  < 2e-16 ***
## s(logDBHcm):SPPSTAGEPIRU 98 1.0000 1.0000  36.534 6.87e-09 ***
## s(logDBHcm):SPPSTAGEPRPE 14 1.0000 1.0000   5.816 0.016800 *  
## s(logDBHcm):SPPSTAGEPRPE 16 1.0000 1.0000   3.494 0.063084 .  
## s(logDBHcm):SPPSTAGEPRPE 17 1.0000 1.0000  23.047 3.09e-06 ***
## s(logDBHcm):SPPSTAGEPRPE 23 1.5500 1.7942   2.852 0.126479    
## s(logDBHcm):SPPSTAGEPRPE 26 1.0000 1.0000  23.426 2.59e-06 ***
## s(logDBHcm):SPPSTAGEPRPE 29 1.0000 1.0000   3.862 0.050810 .  
## s(logHTcm):SPPSTAGEACPE 17  1.0000 1.0000   0.551 0.458696    
## s(logHTcm):SPPSTAGEACPE 26  1.0000 1.0000   0.042 0.838394    
## s(logHTcm):SPPSTAGEACPE 98  1.0000 1.0000   0.019 0.890443    
## s(logHTcm):SPPSTAGEACRU 14  1.0000 1.0000   0.426 0.514697    
## s(logHTcm):SPPSTAGEACRU 16  1.0000 1.0000   0.128 0.720471    
## s(logHTcm):SPPSTAGEACRU 26  1.0000 1.0000   0.437 0.509197    
## s(logHTcm):SPPSTAGEACRU 29  1.0000 1.0000   0.000 0.986487    
## s(logHTcm):SPPSTAGEACSA 14  1.0000 1.0000   0.006 0.938250    
## s(logHTcm):SPPSTAGEACSA 16  0.6146 0.6146  35.906 4.89e-06 ***
## s(logHTcm):SPPSTAGEACSA 17  1.8681 1.9836   3.303 0.052826 .  
## s(logHTcm):SPPSTAGEACSA 26  1.0000 1.0000   0.057 0.812109    
## s(logHTcm):SPPSTAGEACSA 29  1.0000 1.0000   0.295 0.587747    
## s(logHTcm):SPPSTAGEACSA 98  1.0000 1.0000   0.713 0.399588    
## s(logHTcm):SPPSTAGEBEAL 14  1.0000 1.0000   1.106 0.294254    
## s(logHTcm):SPPSTAGEBEAL 16  1.0000 1.0000   0.101 0.750507    
## s(logHTcm):SPPSTAGEBEAL 17  2.8624 2.9821  12.693 4.03e-07 ***
## s(logHTcm):SPPSTAGEBEAL 23  1.0000 1.0000   0.081 0.776357    
## s(logHTcm):SPPSTAGEBEAL 26  1.5280 1.8050   1.201 0.399854    
## s(logHTcm):SPPSTAGEBEAL 29  1.0000 1.0000   0.097 0.755601    
## s(logHTcm):SPPSTAGEBEAL 98  1.5447 1.7934   0.799 0.345325    
## s(logHTcm):SPPSTAGEBEPA 14  1.0000 1.0000   0.125 0.724112    
## s(logHTcm):SPPSTAGEBEPA 16  1.0000 1.0000   2.010 0.157857    
## s(logHTcm):SPPSTAGEBEPA 17  1.1492 1.2780   0.259 0.584128    
## s(logHTcm):SPPSTAGEBEPA 23  1.0000 1.0000   0.065 0.799534    
## s(logHTcm):SPPSTAGEBEPA 26  1.0000 1.0000   0.492 0.483732    
## s(logHTcm):SPPSTAGEBEPA 29  1.0000 1.0000   0.044 0.834086    
## s(logHTcm):SPPSTAGEFAGR 14  1.0000 1.0000   0.606 0.437345    
## s(logHTcm):SPPSTAGEFAGR 16  1.0000 1.0000   6.111 0.014285 *  
## s(logHTcm):SPPSTAGEFAGR 17  1.1686 1.2967   1.387 0.415472    
## s(logHTcm):SPPSTAGEFAGR 26  1.0000 1.0000   0.043 0.835621    
## s(logHTcm):SPPSTAGEFAGR 29  0.6781 0.6781  62.894 4.79e-10 ***
## s(logHTcm):SPPSTAGEFAGR 98  2.0526 2.2926  17.211 2.04e-07 ***
## s(logHTcm):SPPSTAGEPIRU 98  2.5078 2.8081   2.889 0.024406 *  
## s(logHTcm):SPPSTAGEPRPE 14  1.0000 1.0000   0.194 0.660419    
## s(logHTcm):SPPSTAGEPRPE 16  1.0000 1.0000  31.753 5.75e-08 ***
## s(logHTcm):SPPSTAGEPRPE 17  1.0000 1.0000  41.740 7.09e-10 ***
## s(logHTcm):SPPSTAGEPRPE 23  1.0061 1.0088   0.065 0.851093    
## s(logHTcm):SPPSTAGEPRPE 26  1.5285 1.7808   0.313 0.720657    
## s(logHTcm):SPPSTAGEPRPE 29  1.0000 1.0000   0.134 0.715177    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Rank: 271/290
## R-sq.(adj) =  0.994   Deviance explained = 99.6%
## GCV = 0.063204  Scale est. = 0.038352  n = 316

Plotting the model:

plot(allom_model, scale=0)

## doesn't seem to be any discernible pattern, however it is hard to determine with 42 plots

## Plotting predicted values, observed values, and confidence interval for each combination of spp and standage:

Write function to plot predictions vs observed for all combinations of SPP and STANDAGE:

plot_pred <- function(data, model, pred.length = 50, plotcov = c("logDBHcm, logHTcm")) {
  ## Build predicted values into dataframe:
  numspp <- length(levels(data$SPP4))
  numstage <- length(levels(data$STANDAGE))
  
  ## write in STANDAGE levels for each SPP
  stagevec <- rep(levels(data$STANDAGE), numspp)
  
  ## write in SPP
  sppvec <- character()
  
  for (i in 1:numspp) {
    newspp <- rep(levels(data$SPP4)[i], numstage)
    sppvec <- c(sppvec, newspp)
  }
  ## combine into dataframe
  fac_data <- data.frame(sppvec, stagevec)
  
  ## get x scale for each STANDAGE and covariate
    STAGE <- levels(data$STANDAGE)
    minDBH <- numeric(length = length(levels(data$STANDAGE)))
    maxDBH <- numeric(length = length(levels(data$STANDAGE)))
    minHT <- numeric(length = length(levels(data$STANDAGE)))
    maxHT <- numeric(length = length(levels(data$STANDAGE)))
    ranges <- data.frame(STAGE, minDBH, maxDBH, minHT, maxHT)
  
  for (i in 1:length(levels(data$STANDAGE))) {
    newsub <- subset(data, STANDAGE == levels(data$STANDAGE)[i])
    ranges$minDBH[i] <- min(newsub$logDBHcm)
    ranges$maxDBH[i] <- max(newsub$logDBHcm)
    ranges$minHT[i] <- min(newsub$logHTcm)
    ranges$maxHT[i] <- max(newsub$logHTcm)
  }
  
  ## write in regularly spaces covariate values to predict upon
  logDBHcm <- numeric()
  logHTcm <- numeric()
  for (i in 1:nrow(fac_data)) {
    range <- subset(ranges, STAGE == stagevec[i])
    logDBHcm <- c(logDBHcm, seq(range$minDBH, range$maxDBH, length.out = pred.length))
    logHTcm <- c(logHTcm, seq(range$minHT, range$maxHT, length.out = pred.length))
  }
  
  SPP4 <- rep(fac_data$sppvec, each = pred.length)
  STANDAGE <- rep(fac_data$stagevec, each = pred.length)
  
  new_data <- data.frame(SPP4, STANDAGE, logDBHcm, logHTcm)
  
  ## Create SPPSTAGE variable in new data for prediction
  for (i in 1:length(new_data$SPP4)) {
    new_data$SPPSTAGE[i] <- paste(new_data$SPP4[i], new_data$STANDAGE[i])
  }
  new_data$SPPSTAGE <- as.factor(new_data$SPPSTAGE)
  
  ## Eliminate covariate combinations in predictive data for which there was no observed data
  del_vec <- character()
  for (i in 1:nrow(new_data)) {
    if (!any(grepl(new_data$SPPSTAGE[i], levels(data$SPPSTAGE)))) {
      new <- rownames(new_data)[i]
      del_vec <- c(del_vec, new)
    }
    else {}
  }
  
  del_vec <- as.numeric(del_vec[!is.na(del_vec)])
  new_data <- new_data[-del_vec,]
  
  ## write predictions to object
  predicted_data <- predict(model, newdata = new_data, se.fit = TRUE)
  
  ## write in predictions and standard errors to dataframe
  new_data$logABOVEg <- predicted_data$fit
  new_data$SE <- predicted_data$se.fit
  new_data$CI.low <- new_data$logABOVEg - predicted_data$se.fit
  new_data$CI.high <- new_data$logABOVEg + predicted_data$se.fit
  
  plots <- list()
  if (grepl("logDBHcm", plotcov) && grepl("logHTcm", plotcov)) {
    DBHplot <- ggplot(new_data, aes(logDBHcm, logABOVEg)) +
                geom_line(size = 1) +
                geom_line(aes(logDBHcm, CI.high), color = "palevioletred1") +
                geom_line(aes(logDBHcm, CI.low), color = "palevioletred1") +
                geom_jitter(data = data, aes(logDBHcm, logABOVEg), size = 0.75) +
                facet_grid(SPP4 ~ STANDAGE, scales = "free")
    
    HTplot  <- ggplot(new_data, aes(logHTcm, logABOVEg)) +
                geom_line(size = 1) +
                geom_line(aes(logHTcm, CI.high), color = "palevioletred1") +
                geom_line(aes(logHTcm, CI.low), color = "palevioletred1") +
                geom_jitter(data = data, aes(logHTcm, logABOVEg), size = 0.75) +
                facet_grid(SPP4 ~ STANDAGE, scales = "free")
    print(DBHplot)
    print(HTplot)
  }
  else if (grepl("logHTcm", plotcov) && !grepl("logDBHcm", plotcov)) {
    HTplot  <- ggplot(new_data, aes(logHTcm, logABOVEg)) +
                geom_line(size = 1) +
                geom_line(aes(logHTcm, CI.high), color = "palevioletred1") +
                geom_line(aes(logHTcm, CI.low), color = "palevioletred1") +
                geom_jitter(data = allom_data_trunc, aes(logHTcm, logABOVEg), size = 0.75) +
                facet_grid(SPP4 ~ STANDAGE, scales = "free")
    print(HTplot)
  }
  else if (grepl("logDBHcm", plotcov) && !grepl("logHTcm", plotcov)) {
    DBHplot <- ggplot(new_data, aes(logDBHcm, logABOVEg)) +
                geom_line(size = 1) +
                geom_line(aes(logDBHcm, CI.high), color = "palevioletred1") +
                geom_line(aes(logDBHcm, CI.low), color = "palevioletred1") +
                geom_jitter(data = data, aes(logDBHcm, logABOVEg), size = 0.75) +
                facet_grid(SPP4 ~ STANDAGE, scales = "free")
    print(DBHplot)
  }
  else {
    print("error: plot cov must contain either logHTcm, logDBHcm or both")
  }
}

plot_pred(data = allom_data_trunc, model = allom_model)

anova(allom_model)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## logABOVEg ~ s(logDBHcm, by = SPPSTAGE, k = 4) + s(logHTcm, by = SPPSTAGE, 
##     k = 4) + SPP4 * STANDAGE
## 
## Parametric Terms:
##               df     F  p-value
## SPP4           7 0.992    0.438
## STANDAGE       6 2.578    0.020
## SPP4:STANDAGE 30 4.586 3.36e-11
## 
## Approximate significance of smooth terms:
##                                edf Ref.df       F  p-value
## s(logDBHcm):SPPSTAGEACPE 17 1.0000 1.0000  12.602 0.000482
## s(logDBHcm):SPPSTAGEACPE 26 1.0000 1.0000  28.241 2.81e-07
## s(logDBHcm):SPPSTAGEACPE 98 1.0000 1.0000  28.734 2.25e-07
## s(logDBHcm):SPPSTAGEACRU 14 1.0000 1.0000  14.472 0.000190
## s(logDBHcm):SPPSTAGEACRU 16 1.0000 1.0000   6.096 0.014402
## s(logDBHcm):SPPSTAGEACRU 26 1.0000 1.0000   2.735 0.099759
## s(logDBHcm):SPPSTAGEACRU 29 1.0000 1.0000   4.961 0.027058
## s(logDBHcm):SPPSTAGEACSA 14 1.0000 1.0000  18.927 2.17e-05
## s(logDBHcm):SPPSTAGEACSA 16 0.4606 0.4606   0.025 0.915474
## s(logDBHcm):SPPSTAGEACSA 17 1.0000 1.0000  10.772 0.001219
## s(logDBHcm):SPPSTAGEACSA 26 1.0000 1.0000  17.680 3.95e-05
## s(logDBHcm):SPPSTAGEACSA 29 1.0000 1.0000   1.589 0.208941
## s(logDBHcm):SPPSTAGEACSA 98 1.4033 1.6758 155.808  < 2e-16
## s(logDBHcm):SPPSTAGEBEAL 14 1.0000 1.0000   8.167 0.004727
## s(logDBHcm):SPPSTAGEBEAL 16 1.0000 1.0000   0.038 0.845011
## s(logDBHcm):SPPSTAGEBEAL 17 1.0000 1.0000  14.943 0.000150
## s(logDBHcm):SPPSTAGEBEAL 23 1.0000 1.0000   1.524 0.218468
## s(logDBHcm):SPPSTAGEBEAL 26 1.0000 1.0000  13.976 0.000242
## s(logDBHcm):SPPSTAGEBEAL 29 1.0000 1.0000   4.091 0.044488
## s(logDBHcm):SPPSTAGEBEAL 98 1.4009 1.6002 197.387  < 2e-16
## s(logDBHcm):SPPSTAGEBEPA 14 1.0000 1.0000   2.361 0.126018
## s(logDBHcm):SPPSTAGEBEPA 16 1.0000 1.0000   0.255 0.614217
## s(logDBHcm):SPPSTAGEBEPA 17 1.0000 1.0000  50.511 1.67e-11
## s(logDBHcm):SPPSTAGEBEPA 23 1.0000 1.0000  26.291 6.86e-07
## s(logDBHcm):SPPSTAGEBEPA 26 1.0000 1.0000  51.077 1.31e-11
## s(logDBHcm):SPPSTAGEBEPA 29 1.0000 1.0000   0.078 0.780913
## s(logDBHcm):SPPSTAGEFAGR 14 1.0000 1.0000   6.648 0.010661
## s(logDBHcm):SPPSTAGEFAGR 16 1.0000 1.0000   5.996 0.015222
## s(logDBHcm):SPPSTAGEFAGR 17 1.6803 1.8858  23.388 2.36e-07
## s(logDBHcm):SPPSTAGEFAGR 26 1.0000 1.0000  39.876 1.59e-09
## s(logDBHcm):SPPSTAGEFAGR 29 0.3243 0.3243  61.359 1.36e-05
## s(logDBHcm):SPPSTAGEFAGR 98 1.0000 1.0000 361.507  < 2e-16
## s(logDBHcm):SPPSTAGEPIRU 98 1.0000 1.0000  36.534 6.87e-09
## s(logDBHcm):SPPSTAGEPRPE 14 1.0000 1.0000   5.816 0.016800
## s(logDBHcm):SPPSTAGEPRPE 16 1.0000 1.0000   3.494 0.063084
## s(logDBHcm):SPPSTAGEPRPE 17 1.0000 1.0000  23.047 3.09e-06
## s(logDBHcm):SPPSTAGEPRPE 23 1.5500 1.7942   2.852 0.126479
## s(logDBHcm):SPPSTAGEPRPE 26 1.0000 1.0000  23.426 2.59e-06
## s(logDBHcm):SPPSTAGEPRPE 29 1.0000 1.0000   3.862 0.050810
## s(logHTcm):SPPSTAGEACPE 17  1.0000 1.0000   0.551 0.458696
## s(logHTcm):SPPSTAGEACPE 26  1.0000 1.0000   0.042 0.838394
## s(logHTcm):SPPSTAGEACPE 98  1.0000 1.0000   0.019 0.890443
## s(logHTcm):SPPSTAGEACRU 14  1.0000 1.0000   0.426 0.514697
## s(logHTcm):SPPSTAGEACRU 16  1.0000 1.0000   0.128 0.720471
## s(logHTcm):SPPSTAGEACRU 26  1.0000 1.0000   0.437 0.509197
## s(logHTcm):SPPSTAGEACRU 29  1.0000 1.0000   0.000 0.986487
## s(logHTcm):SPPSTAGEACSA 14  1.0000 1.0000   0.006 0.938250
## s(logHTcm):SPPSTAGEACSA 16  0.6146 0.6146  35.906 4.89e-06
## s(logHTcm):SPPSTAGEACSA 17  1.8681 1.9836   3.303 0.052826
## s(logHTcm):SPPSTAGEACSA 26  1.0000 1.0000   0.057 0.812109
## s(logHTcm):SPPSTAGEACSA 29  1.0000 1.0000   0.295 0.587747
## s(logHTcm):SPPSTAGEACSA 98  1.0000 1.0000   0.713 0.399588
## s(logHTcm):SPPSTAGEBEAL 14  1.0000 1.0000   1.106 0.294254
## s(logHTcm):SPPSTAGEBEAL 16  1.0000 1.0000   0.101 0.750507
## s(logHTcm):SPPSTAGEBEAL 17  2.8624 2.9821  12.693 4.03e-07
## s(logHTcm):SPPSTAGEBEAL 23  1.0000 1.0000   0.081 0.776357
## s(logHTcm):SPPSTAGEBEAL 26  1.5280 1.8050   1.201 0.399854
## s(logHTcm):SPPSTAGEBEAL 29  1.0000 1.0000   0.097 0.755601
## s(logHTcm):SPPSTAGEBEAL 98  1.5447 1.7934   0.799 0.345325
## s(logHTcm):SPPSTAGEBEPA 14  1.0000 1.0000   0.125 0.724112
## s(logHTcm):SPPSTAGEBEPA 16  1.0000 1.0000   2.010 0.157857
## s(logHTcm):SPPSTAGEBEPA 17  1.1492 1.2780   0.259 0.584128
## s(logHTcm):SPPSTAGEBEPA 23  1.0000 1.0000   0.065 0.799534
## s(logHTcm):SPPSTAGEBEPA 26  1.0000 1.0000   0.492 0.483732
## s(logHTcm):SPPSTAGEBEPA 29  1.0000 1.0000   0.044 0.834086
## s(logHTcm):SPPSTAGEFAGR 14  1.0000 1.0000   0.606 0.437345
## s(logHTcm):SPPSTAGEFAGR 16  1.0000 1.0000   6.111 0.014285
## s(logHTcm):SPPSTAGEFAGR 17  1.1686 1.2967   1.387 0.415472
## s(logHTcm):SPPSTAGEFAGR 26  1.0000 1.0000   0.043 0.835621
## s(logHTcm):SPPSTAGEFAGR 29  0.6781 0.6781  62.894 4.79e-10
## s(logHTcm):SPPSTAGEFAGR 98  2.0526 2.2926  17.211 2.04e-07
## s(logHTcm):SPPSTAGEPIRU 98  2.5078 2.8081   2.889 0.024406
## s(logHTcm):SPPSTAGEPRPE 14  1.0000 1.0000   0.194 0.660419
## s(logHTcm):SPPSTAGEPRPE 16  1.0000 1.0000  31.753 5.75e-08
## s(logHTcm):SPPSTAGEPRPE 17  1.0000 1.0000  41.740 7.09e-10
## s(logHTcm):SPPSTAGEPRPE 23  1.0061 1.0088   0.065 0.851093
## s(logHTcm):SPPSTAGEPRPE 26  1.5285 1.7808   0.313 0.720657
## s(logHTcm):SPPSTAGEPRPE 29  1.0000 1.0000   0.134 0.715177
## looks like the SPP4 parametric terms arent significant.  STANDAGE and interactions are both significant.

allom_p.coeff <- data.frame(summary(allom_model)$p.coeff)

coefnames <- as.factor(attributes(summary(allom_model)$p.coeff)$names)
coefvals <- summary(allom_model)$p.coeff
attributes(coefvals) <- NULL

allom_coef <- data.frame(coefnames, coefvals)

allom_coefSTANDAGE <- allom_coef[9:14, ]
allom_coefSPP4 <- allom_coef[2:8, ]

## there appears to be some sort of pattern here
ggplot(allom_coefSTANDAGE, aes(x = coefnames, y = coefvals)) +
  geom_point()  

ggplot(allom_coefSPP4, aes(x = coefnames, y = coefvals)) + 
  geom_point()

## x axis is illegible, but could be useful for seeing larger trends
ggplot(allom_coef, aes(x = coefnames, y = coefvals)) + 
  geom_point()

with STANDAGE as class numeric (IGNORE FOR NOW)

# allom_p.coeff <- data.frame(summary(allom_model)$p.coeff)
# 
# coefnames <- as.factor(attributes(summary(allom_model)$p.coeff)$names)
# coefvals <- summary(allom_model)$p.coeff
# attributes(coefvals) <- NULL
# 
# allom_coef <- data.frame(coefnames, coefvals)
# 
# allom_coefSPP4 <- allom_coef[2:8, ]
# allom_coefSTANDAGE <- allom_coef[9:14, ]
# allom_coefSTANDAGE$coefnames <- as.character(allom_coefSTANDAGE$coefnames)
# 
# for (i in 1:nrow(allom_coefSTANDAGE)) {
#   allom_coefSTANDAGE[i, 1] <- gsub("STANDAGE", "", allom_coefSTANDAGE[i, 1])
# }
# 
# allom_coefSTANDAGE$coefnames <- as.numeric(allom_coefSTANDAGE$coefnames)
# allom_coefSTANDAGE
# 
# ggplot(allom_coefSTANDAGE, aes(x = coefnames, y = coefvals)) +
#   geom_point() 
# 
# ggplot(allom_coefSPP4, aes(x = coefnames, y = coefvals)) + 
#   geom_point()
# 
# ## x axis is illegible, but could be useful for seeing larger trends
# ggplot(allom_coef, aes(x = coefnames, y = coefvals)) + 
#   geom_point()